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Abstract 

A  previous  paper  by  the  authors  [FH93]  noted  that  there  was  a  strong  tendency  to 
obtain  near-repeated  knots  in  their  algorithm  for  least  squares  approximation  of  scat- 
tered data  by  multiquadric  functions.  In  this  paper  we  observe  that  this  leads  naturally 
to  the  inclusion  of  derivatives  of  the  multiquadric  basis  function  in  the  approximation, 
and  give  an  algorithm  for  accomplishing  this.  A  comparison  of  the  results  obtained 
with  this  algorithm  and  the  previous  one  is  made.  While  the  multiple  knot  algorithm 
usually  has  the  advantage  in  terms  of  accuracy  and  computational  stability,  there  are 
datasets  for  which  this  is  reversed. 
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1      Introduction 

This  paper  continues  an  investigation  into  the  use  of  multiquadric  functions  to  approximate 
scattered  data.  An  initial  report  on  our  efforts  appears  in  [FH93],  where  we  reported  that 
a  phenomenon  that  seemed  particularly  interesting  was  the  tendency  of  the  optimization 
process  to  result  in  near-repeated  knots.  The  implication  is  that  the  method  was  attempting 
to  build  a  directional  derivative  of  the  basis  function  into  the  approximation.  This  paper 
documents  our  experiences  with  an  algorithm  that  detects  the  occurrence  of  near-repeated 
knots,  and  when  they  occur,  replaces  basis  functions  corresponding  to  near-repeated  knots 
with  a  single  multiquadric  basis  function  plus  appropriate  derivatives  of  the  multiquadric 
basis  function,  at  the  pertinent  knot. 

While  it  is  useful  to  have  the  previous  paper  at  hand,  for  completeness  we  will  briefly 
review  the  background  necessary  to  make  this  report  somewhat  self-contained.  We  restrict 
our  discussion  to  functions  of  two  independent  variables,  the  methods  are  easily  extendible 
to  arbitrary  dimensions,  and  we  expect  that  many  of  the  conclusions  will  carry  over. 

The  scattered  data  approximation  problem  is  easily  described  and  occurs  frequently  in 
many  branches  of  science.  The  problem  occurs  in  any  discipline  where  measurements  are 
taken  at  irregularly  spaced  values  of  two  or  more  independent  variables,  and  is  especially 
prevalent  in  environmental  sciences.  We  will  suppose  that  triples  of  data,(xj,  y},  Zj),  j  = 
1,---,7V  are  given,  assumed  to  be  measurements  (perhaps  with  error)  of  an  underlying 
function  z  =  f(x,y).  The  function  /  is  to  be  approximated  by  a  function  F(x,y)  from  the 
given  data.  A  recent  survey  of  such  methods  is  given  in  [FN91]. 

Multiquadric  functions  were  introduced  for  interpolation  of  scattered  data  by  Hardy 
[HA71];  also  see  [HA90]  for  a  historical  survey  and  many  references.  The  method  is  one 
of  a  class  of  methods  known  now  as  "radial  basis  function  methods"  that  includes  other 
attractive  schemes  such  as  thin  plate  splines  [HD72,  DU76,  and  others].  The  basic  idea 
of  such  methods  is  quite  simple,  and  we  describe  it  in  some  generality;  for  purposes  of 
being  definite  it  is  pertinent  to  note  that  for  the  multiquadric  method  the  radial  function 
is  h(d)  =  J{d2  +  r2).  In  general,  suppose  a  function  of  one  variable,  /?((/),  where  d  denotes 
distance,  is  given. 

For  interpolation  (that  is,  exact  matching  of  the  given  data),  a  basis  function,  Bj(x,y)  = 
h(dj)  is  associated  with  each  data  point.  Here  d-,  =  \/{{x  —  x3)2  +  (y  —  yj)2),  the  distance 
from  (x,y)  to  (x^y^).  Thus  each  basis  function  is  a  translate  of  the  radial  function,  h.  The 
approximation  is  a  linear  combination  of  the  basis  functions,  along  with  some  polynomial 
terms  that  may  be  necessary  in  some  cases,  or  may  be  used  to  assure  that  the  approximation 
method  has  polynomial  precision.  Thus, 

N  M 

F(x,y)  =  Y,<*jBj{x,y)  +  ^6j9j(.r,y)  (1) 

where  {qj}  is  a  set  of  M  polynomials  forming  a  basis  for  polynomials  of  degree  <  m.  The  co- 
efficients a-j  and  bj  are  determined  by  the  linear  system  of  equations  prescribing  interpolation 


of  the  data,  and  exactness  for  polynomials  of  degree  <  m: 

N  M 

J2aJBi(xi>yi)  +Ylbjqj(xi,yi)  =  zt,  i  -  1,---,A' 

*J>  i""  (2) 

Yl ai<li(xJi Vi)  =  °>  «  =  1, ■  •  •  ,M. 

For  multiquadric  basis  functions,  this  system  of  equations  is  known  to  have  a  unique  solution 
for  distinct  (x^, y3)  data  (see,  for  example,  [MI86]);  while  m  may  be  taken  as  zero  (no 
polynomial  terms),  Micchelli's  results  show  the  multiquadric  basis  function  is  positive  definite 
of  order  one,  and  thus  a  constant  term  should  be  included.  We  have  done  so  in  all  our  work.  If 
higher  degree  polynomial  precision  is  desired,  inclusion  of  those  terms  imposes  no  particular 
burden. 

While  interpolation  theory  is  important  and  indicates  something  about  the  suitability 
of  the  class  of  functions  for  approximation  purposes,  our  emphasis  here  is  on  least  squares 
approximation.  This  implies  using  fewer  basis  functions  than  there  are  data  points.  In 
analogy  with  univariate  cubic  splines,  it  is  convenient  to  refer  to  the  points  at  which  the 
radial  basis  functions  are  centered  as  "knots1',  as  was  done  in  [MF92],  and  we  do  so  here. 
If  a  set  of  knot  points,  (uk,Vk),k  =  1,--,A',  with  K  <  N  have  been  specified,  then  the 
problem  of  fitting  a  multiquadric  function  by  least  squares  is  similar  to  that  of  solving 
the  system  of  equations  corresponding  to  those  above  in  the  least  squares  sense.  We  give 
the  details.  Now,  let  Bk{x,y)  denote  the  radial  basis  function  associated  with  the  point 
(iik,Vk),  Bk(x,  y)  =  \/({x  —  Uk)2  +  (y  —  vk)2  +  r2).  The  system  of  equations,  specialized  for 
our  case,  is  now  of  the  form 

K 

^akBk(xt,yt)  +  c  =  z,,  i=  l,--,N 

¥  (3) 

There  is  a  question  of  how  to  treat  the  last  equation,  which  guarantees  precision  for  con- 
stants. In  [MF92]  the  corresponding  constraint  equations  were  imposed  exactly,  rather  than 
approximately,  because  of  physical  considerations.  While  there  is  not  the  corresponding 
physical  situation  here,  we  have  also  imposed  the  last  equation  as  a  constraint.  This  con- 
straint can  be  used  to  reduce  the  size  of  the  system  by  solving  for  a^-  in  terms  of  the  other 
ak  and  substituting  into  the  first  set  of  equations. 

If  the  knot  points  are  a  subset  of  the  data  points,  then  the  system  of  equations  (by 
virtue  of  containing,  as  a  subset,  interpolation  equations)  has  full  rank,  and  thus  guarantees 
a  unique  solution  of  the  least  squares  problem.  When  the  knot  locations  may  differ  from  the 
data  points,  the  problem  of  whether  the  coefficient  matrix  is  of  full  rank  or  not  is  unknown  to 
us,  although  we  feel  certain  that  the  matrix  is  of  full  rank  when  the  knot  points  are  distinct, 
and  have  encountered  no  situations  that  indicate  otherwise. 

The  impetus  behind  our  original  investigation  was  to  obtain  surface  approximations  that 
would  be  efficient  in  subsequent  applications.    That  is,  we  consider  it  to  be  acceptable  to 


expend  considerable  computational  resources  to  obtain  the  approximation  in  a  preprocessing 
step.  Once  obtained,  the  approximation  can  be  evaluated  efficiently,  so  it  would  be  feasible 
to  use  it  numerous  times  in  an  application  program. 

In  the  previous  paper  [FH93]  we  found  that  there  was  a  tendency  of  the  algorithm  to 
converge  on  solutions  where  the  knot  points  were  often  close  to  each  other,  or  near-repeated. 
This  occurred  both  in  clusters  of  two  and  three  knots.  At  that  time,  no  attempt  was  made  to 
made  to  investigate  this  behavior  further.  The  outstanding  result  of  that  reference  was  the 
finding  that  allowing  the  knots  (and  r  value)  to  be  determined  as  part  of  the  minimization 
process  yielded  approximations  that  were  much  more  efficient  than  those  obtained  either  by 
determining  the  knots  a  priori  or  adaptively  (we  outlined  a  greedy  algorithm  for  sequentially 
determining  knot  locations  that  worked  well).  While  we  found  that  using  variable  r  values 
also  yielded  good  improvement,  we  do  not  pursue  that  idea  here. 

Because  the  greedy  algorithm  we  developed  previously  was  generally  used  as  the  starting 
point  for  the  work  reported  here,  we  give  a  brief  outline  of  it. 

a)  Obtain  the  least  squares  fit  by  a  constant  function,  the  average  of  the  data  values.  The 

two  data  points  having  maximum  positive  and  maximum  negative  error  are  taken  to 
be  the  first  two  knots,  (ui,t>i)  and  (u2,t'2)-  The  knot  counter  A  is  set  to  2. 

b)  The  least  squares  multiquadric  fit  with  A'  knots  is  obtained,  and  the  residuals  are 

computed. 

c)  The  maximum  absolute  value  of  the  residuals  is  found  and  the  location  of  this  residual, 

subject  to  the  minimum  knot  separation  value,  is  taken  to  be  the  next  knot  location 
{uKivk)-  At  this  point  the  algorithm  proceeds  to  step  b  unless  the  maximum  number 
of  knot  locations  to  be  computed  has  been  reached. 

As  previously,  we  have  used  a  QRP1  decomposition  of  the  coefficient  matrix  to  solve  the 
least  squares  problem  for  given  knots  and  r  value.  This  provides  a  stable  and  efficient  means 
for  solution  of  the  problem  with  an  indication  if  a  matrix  of  less  than  full  rank  is  encountered. 

In  order  to  test  the  algorithms  we  have  used  a  number  of  data  sets.  Several  of  these  are 
based  on  previously  published  and  widely  available  {x,y)  data  sets  and  parent  functions.  We 
have  also  used  a  few  less  readily  available  data  sets  that  we  are  willing  to  share  with  anyone 
interested  in  obtaining  them.  Table  1  gives  a  summary  of  most  of  the  data  sets. 

As  mentioned  in  the  opening  paragraph,  the  primary  purpose  of  this  paper  is  to  discuss 
some  further  results  we  have  obtained  in  allowing  the  coalescence  of  knots.  A  discussion  of 
the  algorithm  we  used  is  given  in  Section  2.  Results  of  calculations  obtained  using  several 
data  sets  and  comparison  with  the  previous  method  where  only  simple  (but  possibly  near- 
repeated)  knots  are  allowed,  is  given  in  Section  3.  Section  4  summarizes  our  experiences. 


n.m  This  refers  to  point  set  n  and  function  ???  from  [FR82],  for  n  =1,  2,  and  3,  and 
m  =  l,---,6.  n  =  1  is  25  points,  n  =  2  is  33  points,  n  =  3  is  100  points,  n  =  4 
refers  to  the  200  point  data  set  used  in  [MF92].  m=l  is  the  humps  and  dip  function, 
m  =  2  is  the  cliff,  m  =  3  is  the  saddle,  m  =  4  is  the  gentle  hill,  m  =  5  is  the  steep  hill, 
and  m  =  6  is  the  sphere.  In  addition,  m  =  7  refers  to  the  curved  valley  function  from 
[NI78]. 

GT  This  refers  to  the  thinned  glacier  data  consisting  of  678  points,  with  certain  contour 
lines  removed,  from  [MF92]. 

GL  This  refers  to  the  thinned  glacier  data  consisting  of  873  points. 

HF  This  is  the  data  set  from  [MF92]  generated  to  have  point  density  approximately  pro- 
portional to  curvature,  consisting  of  500  points. 


Table  1:  Data  Sets  Used  Extensively  in  Tests 

2      Derivatives   of  Multiquadric  Basis   Functions   and 
Multiple  Knots 

Once  again  we  have  used  Matlab  '  and  LEASTSQ  (a  Levinburg-Marquardt  method)  to  solve 
the  nonlinear  minimization  problem 

N  K 

min  J^fe  -  51  akBk{xi,yi)  -  c)2  (4) 

where  the  minimization  is  taken  over  all  (u*,t'fc),  r,  the  a^,  and  c  (with  the  last  equation  of 
(3)  imposed  as  a  constraint)  as  an  equivalent  two  step  minimization.  Thus,  for  each  given 
knot  configuration  and  r  value,  the  least  squares  solution  of  (3)  is  computed  as  a  step  toward 
(4).  This  results  in  the  solution  of  a  simpler,  but  equivalent  problem  since  2K  parameters 
are  eliminated  from  (4)  by  imposing  the  condition  that  the  values  of  the  a,j  and  c  be  always 
taken  as  obtained  from  the  least  squares  solution  of  (3).  Hence,  our  final  process  is  more 
properly  written  as 

N  N 

minminJ^Zi  -  £^  akBk{xi,  y,-)  -  c]2  (5) 

1=1  A:=l 

where  the  inner  minimization  is  over  the  a,j  and  c  (least  squares  solution  of  (3)),  and  the 
outer  minimization  is  over  the  knot  locations  and  the  value  of  the  parameter  r.  The  global 
minimum  of  each  of  the  two  problems  are  clearly  the  same.  Eq.  (5)  is  the  more  restrictive, 
but  any  minimum  of  (4)  is  a  local  minimum  of  (5),  else  a  better  solution  is  attainable  for  (4). 


1  Math  Works,  24  Prime  Park  Way,  Natick.  MA  01760 


This  does  not  imply  that  the  iterative  methods  employed  to  solve  (5)  would  work  equally 
well,  nor  find  the  same  local  minima,  when  applied  to  (4). 

When  two  knots  approach  each  other,  the  system  (3)  becomes  ill-  conditioned,  and  the 
coefficients  for  the  two  knots  tend  to  become  large  and  of  opposite  sign.  This  reminds  one  of 
derivatives,  in  this  case  a  directional  derivative.  Because  the  two  coefficients  are  not  of  the 
same  magnitude,  the  appropriate  replacement  basis  functions  seem  to  be  one  of  the  basis 
functions,  and  a  directional  derivative  of  the  basis  function.  In  practice,  the  directional 
derivative  is  represented  as  a  linear  combination  of  the  partial  derivatives  with  respect  to  x 
and  to  y.  To  prove  this  is  the  proper  change  requires  showing  that  the  coefficients  grow  at  no 
greater  a  rate  than  the  reciprocal  of  the  distance  between  the  knots.  The  work  of  Narcowich 
and  Ward  [NW91]  shows  that  the  norm  of  the  inverse  of  the  multiquadric  interpolation  matrix 
grows  no  more  rapidly  than  linear  with  the  reciprocal  of  the  minimum  distance  between  data 
points.  It  seems  likely  a  similar  result  holds  for  the  norm  of  the  pseudoinverse  of  the  system 
matrix  for  (3).  We  have  verified  this  computationally  for  some  cases. 

The  contribution  of  the  two  terms  corresponding  to  the  coalescent  knots  (say  the  jth 
and  the  kth)  are  of  the  form  ajBj(x,y)  +  akBk(x,y).  The  six  parameters,  Oj,  a^,  i/j,  Vj, 
Uk,  and  Vf.,  in  those  terms  are  replaced  by  the  five  parameters  in  the  terms  a3Bj(x,y)  -f 
b3(Bj(x,y))x  +  Cj(Bj(x,y))y.  In  the  previous  expression  the  subscripts  x  and  y  refer  to 
partial  differentiation.  We  will  refer  to  such  a  knot  as  a  double  knot. 

Now,  as  it  turns  out,  the  coalescence  of  two  knots,  while  the  most  prevalent,  was  by  no 
means  the  only  situation  we  encountered.  The  coalescence  of  three  knots  (or,  once  two  knots 
have  coalesced,  the  approach  of  a  third  knot  to  the  location  of  a  double  knot)  then  results  in 
the  occurrence  of  a  second  derivative,  a  directional  derivative  of  a  directional  derivative  (the 
directions  not  necessarily  being  the  same).  In  practice  the  replacement  of  the  basis  functions 
for  the  ith,  jth,  and  kth  knots  is  by  the  linear  combination  a^Bj  +  bj(Bj)x  -f  Cj{B3)y  -f 
dj(Bj)xx  +  ej(Bj)xy  -f  fj{Bj)yy.  While  we  have  instances  of  near-coalescence  of  four  knots, 
our  present  code  does  not  attempt  to  handle  such  cases  properly. 

The  algorithm  is  required  to  change  basis  functions  when  it  appears  that  knots  are 
coalescing.  We  have  used  a  simple  idea  that  seems  to  be  effective.  We  iterate  the  solution  of 
(5)  in  stages  having  increasingly  stringent  convergence  tolerances.  As  each  tolerance  level  is 
met,  a  check  for  near-repeated  knots  is  made,  "near"  being  a  changing  tolerance  consistent 
with  the  tolerance  for  the  optimization.  When  pairs  of  knots  are  found  within  the  tolerance, 
the  new  basis  is  generated,  and  the  optimization  process  continued.  When  no  new  repeated 
knots  are  obtained,  the  tolerance  is  decreased  by  a  factor  of  two.  We  have  observed  that  the 
inclusion  of  the  derivative  terms,  vice  the  basis  functions  at  near-  repeated  knots,  generally 
yields  essentially  the  same  approximation,  as  it  should,  for  both  double  and  triple  knots.  The 
sequence  of  tolerance  levels  we  used  is:  0.02,  0.01,  0.005  for  problems  posed  on  the  [0,  l]2 
square. 

3      Approximations  with  Repeated  Knots 

The  results  of  our  investigation  concern  the  ability  of  the  algorithm  to  obtain  good  ap- 
proximations in  comparison  with  the  algorithm  permitting  only  simple  (but  perhaps  near- 


repeated)  knots,  the  degree  of  dependence  on  the  initial  guess  for  three  different  datasets 
(two  in  some  detail),  and  a  comparison  of  condition  numbers  of  the  coefficient  matrices  for 
the  simple  and  multiple  knot  problems. 


dset         rs  rm  rmss       rmsm      grmss     grmsm      knotss         knots, 


3.1 

0.0976 

0.1588 

0.0081 

0.0059 

0.0115 

0.0107 

2nd,lnt 

Id, It, lnt 

3.2 

7e-6 

0.0002 

0.0074 

0.0073 

0.0186 

0.2477 

3nd,lnt 

3d,lt,lnq 

3.3 

0.3116 

0.2502 

0.0014 

0.0008 

0.0016 

0.0011 

lnt 

It 

3.4 

0.2885 

0.3118 

0.0005 

0.0005 

0.0006 

0.0005 

2nd 

3d 

3.5 

0.3720 

0.3886 

0.0008 

0.0007 

0.0010 

0.0010 

3nd 

2d 

3.6 

2.3280 

2.3044 

0.0005 

0.0005 

0.0007 

0.0007 

3.7 

0.0990 

0.3304 

0.0403 

0.0279 

0.0733 

0.0324 

4nd 

4d 

4.1 

0.0958 

0.1588 

0.0070 

0.0017 

0.0071 

0.0024 

lnt 

It 

4.2 

2e-5 

2e-5 

0.0035 

0.0035 

0.0083 

0.0082 

4nd 

3d,  lnd 

4.3 

0.2577 

0.2513 

0.0002 

0.0002 

0.0002 

0.0002 

4nd 

Id, lnd 

4.4 

0.4163 

0.4625 

7e-6 

6e-6 

7e-6 

6e-6 

lnd 

3d 

4.5 

0.3234 

0.3151 

0.0002 

0.0003 

0.0002 

0.0003 

4nd 

2nd 

4.6 

0.5668 

0.5680 

5e-5 

4e-5 

5e-5 

4e-5 

5nd 

2d,2nt,2nd 

4.7      0.1111     0.1134     0.0115     0.0114     0.0202     0.0287     3nd,lnt  5d,lt 


Table  2:  Comparison  of  Results:  Near-Repeated  Knots  vs.  Multiple  Knots 

We  first  describe  the  layout  of  the  Table  2,  which  gives  a  comparison  of  the  results  of 
the  simple  knot  and  multiple  knot  algorithms  for  fourteen  datasets.  The  results  are  for  the 
surface  fits  obtained  by  starting  from  the  initial  guess  obtained  by  the  greedy  algorithm 
[FH93],  with  r=  0.3,and  closeness  tolerance,  ctol=0.1.  The  3.n  datasets  were  fit  using  twelve 
knot  locations,  while  the  4.n  datasets  were  fit  using  twenty  knot  locations.  The  columns 
of  the  table  contain,  in  pairs,  the  values  of  r,  the  rms  errors  for  the  data  points,  the  rms 
errors  on  a  21x21  grid,  and  information  about  the  near-repeated  or  multiple  knots  obtained. 
The  subscripts  s  and  m  denote  the  simple  knot  and  multiple  knot  algorithms,  respectively. 
The  terms  2d  and  It,  for  example,  refer  to  2  double  knots,  and  1  triple  knot.  Likewise,  lnd 
refers  to  1  near-double  knot.  In  some  cases  two  double  knots  may  have  been  obtained,  with 
another  simple  knot  nearby  one  of  the  double  knots  -  such  a  case  might  be  denoted  by  2d, 
Int. 

For  a  given  dataset,  the  values  of  r  obtained  are  mostly  comparable,  with  larger  values 
generally  being  obtained  by  the  multiple  knot  algorithm.  In  the  instances  where  the  value 
of  r  is  not  larger,  the  values  are  essentially  the  same,  the  exception  being  dataset  3.3.  The 
cases  showing  the  most  difficulty  are  the  cases  3.2  and  4.2  (the  cliff  function),  which  has  a 


tendency  toward  small  values  of  the  parameter  r.  The  effects  of  this  tendency  are  discussed 
in  a  later  paragraph. 

The  rms  errors  are  mostly  comparable,  with  the  errors  obtained  by  the  multiple  knot 
algorithm  generally  smaller,  with  some  ties,  and  one  instance  where  the  errors  are  bigger. 
This  occurs  for  dataset  4.5,  but  the  errors  for  both  approximations  are  quite  small  for  this 
case. 

The  rms  errors  on  the  grid  (grms)  generally  follow  the  same  pattern,  with  one  major 
exception.  On  the  dataset  3.2,  the  errors  for  the  grid  are  quite  large.  Closer  inspection  of 
this  particular  case  reveals  there  is  one  grid  point  (near  the  triple,  almost  quadruple,  knot) 
where  the  surface  is  badly  ill-behaved.  Coefficients  of  the  first  derivative  terms  are  on  the 
order  of  5  to  20,  while  the  coefficients  of  the  two  multiquadric  basis  functions  are  reasonably 
large,  on  the  order  of  103.  However,  the  value  of  r  is  very  small  in  this  case,  about  1.5e-4, 
thus  the  derivatives  of  the  basis  functions  very  rapidly  approach  the  value  one  in  one  of  the 
coordinate  directions  (the  value  is  about  0.99  at  a  distance  of  0.001  from  the  knot).  Hence, 
very  rapid  changes  in  the  value  of  the  function  occur.  While  otherwise  similar  situations 
occur  in  other  cases,  the  key  difference  here  is  the  small  value  of  r  that  occurs. 

In  order  to  discover  whether  the  initial  guess  was  contributing  to  the  problem,  or  whether 
the  tendency  toward  a  small  r  value  and  multiple  knots  was  driven  by  the  data,  we  ran  a 
number  initial  guesses  for  the  dataset  3.2.  The  initial  conditions  were  obtained  from  the 
greedy  algorithm  with  various  values  of  r  and  closeness  tolerance  (ctol).  The  results  of  the 
computations  are  summarized  in  Table  3.  The  table  gives  the  parameter  ctol,  the  initial 
value  of  r,  the  final  value  of  r,  the  rms  errors  at  the  data  points  and  on  the  grid,  and  the 
number,  multiplicity,  and  approximate  location  of  the  multiple  knots.  From  this  we  see  that 
while  there  is  a  tendency  toward  multiple  knots,  no  particular  locations  are  tended  toward, 
although  it  seems  natural  that  repeated  knots  invariably  occur  near  the  diagonal.  The  fact 
that  the  surface  is  a  function  of  one  variable  (that  is,  y-x)  with  rapid  changes  occurring 
near  the  diagonal  may  account  for  this  circumstance.  The  fact  that  the  data  is  not  regular 
may  account  for  a  large  number  of  comparable  local  minima.  One  case  is  of  particular  note: 
The  case  in  the  second  line  involves  three  triple  and  one  double  knot,  however  two  of  the 
triple  knots  are  close  together,  yielding  a  near-repeated  knot  of  multiplicity  six.  This  could 
be  expected  to  lead  to  poor  behavior  of  the  surface.  Closer  inspection  revealed  that,  even 
though  the  coefficients  are  rather  large  (order  107),  the  surface  is  well  behaved  and  yielded 
one  of  the  better  approximations  among  this  group.  This  good  behavior  is  undoubtedly 
enhanced  by  the  fact  that  the  r  value  is  large  for  this  surface.  In  the  three  cases  where 
the  r  value  is  smaller  than  10~4  the  resulting  surface  is  poorly  behaved  near  a  multiple  or 
near-multiple  knot,  this  being  reflected  by  the  larger  grms0  values  shown. 

We  ran  a  number  of  different  initial  guesses  for  case  3.1  to  determine  both  the  robustness 
of  the  optimization  routine  and  to  confirm  the  tendency  of  the  method  to  converge  on 
multiple  knots  from  a  variety  of  initial  guesses.  Nine  different  initial  guesses  were  made, 
six  resulting  from  the  greedy  algorithm  with  different  closeness  tolerances  and  r  values,  and 
three  with  random  initial  knot  locations.  As  is  seen  in  Table  4,  the  results  for  the  various 
initial  guesses  are  strikingly  similar:  except  for  two  cases  the  r  values  are  all  close  to  0.16,  all 
but  one  result  in  a  triple  or  near-triple  knot  near  the  dip  at  about  (.45, .78),  and  all  but  one 
have  comparable  accuracies  for  rms0  and  grms0.  The  one  odd  case  had  the  starting  value 


of  0.01  for  r.   A  wide-ranging  difference  is  that  different  double  and  near-double  knots  may 
occur  for  the  various  runs.  Unlike  dataset  3.2,  this  case  appears  to  be  quite  stable. 


ctol        r  r0  rms0      grms0     multiple  knots 


t@(.09,.14),  3d@(.38,.38,  (.60, .60),  (.89, .88)) 
3t@(.31,.36),  (.76,.79),  (.77,-78),  ld@(.12,.13) 
2d@(.38,.40),  (.67,.66) 
3d@(.14,.14),(.44,.44),  (.68,-67) 
4d@(.13,.ll),(.37,.39),  (.69,-65),  (.92,-91) 
d@(.69,.69),  2nd@(.14,.12),  (.33,-40) 
d@(.44,.44) 
d@(.69,.69),  2nd@(.14,.ll),  (.33,-40) 


0.1 

0.3 

0.0002 

0.0073 

0.2496 

0.25 

0.1 

0.0239 

0.0060 

0.0085 

0.1, 

0.1 

0.0247 

0.0071 

0.0111 

0.2 

0.1 

0.0247 

0.0071 

0.0113 

0.25 

0.01 

0.0739 

0.0054 

0.0095 

0.1 

0.01 

5e-5 

0.0071 

0.0259 

0.2 

0.01 

0.0001 

0.0071 

0.0114 

0.0 

0.01 

5e-5 

0.0071 

0.0260 

Table  3:  Results  with  Different  Initial  Guesses  for  Case  3.2 


ctol 

r 

sk 

r0 

rms0 

grms0 

multiple  knots 

0.0 

0.3 

g 

0.1598 

0.0067 

0.0113 

t@(.43,.76),  2nd@(.86,.23),  (.83, 

.42) 

0.1 

0.3 

g 

0.1588 

0.0059 

0.0107 

t@(.43,.75),  d@(.78,.43),  d  is  nt 

0.2 

0.3 

g 

0.1712 

0.0063 

0.0107 

t@(.44,.74),  d@(.15,.26) 

0.25 

0.3 

g 

0.2013 

0.0075 

0.0107 

t@(.43,.79),  d@(.20,.33) 

0.1 

0.1 

g 

0.1450 

0.0066 

0.0111 

t@(.43,.77) 

0.1 

0.01 

g 

0.0737 

0.0118 

0.0157 

d@(.72,.36) 

- 

0.3 

r 

0.1454 

0.0066 

0.0112 

t@(.43,.77) 

- 

0.3 

r 

0.1566 

0.0061 

0.0077 

t@(.42,.75) 

- 

0.3 

r 

0.1600 

0.0067 

0.0113 

t@(.43,.76),  2nd@(.83,.43),  (.86, 

.23) 

0.3 

r 

0.1527 

0.0067 

0.0111 

nt@(.44,.75) 

Table  4:  Results  with  Different  Initial  Guesses  for  Case  3.1 

We  also  ran  some  different  initial  conditions  on  dataset  3.7.  The  results  in  this  case  seem 
to  indicate  complications  somewhere  between  the  datasets  3.1  and  3.2.  Dataset  3.7  has  a 
tendency  toward  many  multiple  knots,  perhaps  because  it  has  a  more  complex  behavior  of 
the  surface  than  any  of  the  other  datasets. 


The  occurrence  of  near-multiple  knots,  as  we  noted  previously,  must  have  a  unfavorable 
effect  on  the  condition  number  of  the  coefficient  matrix,  and  in  order  to  investigate  this  we 
computed  the  condition  number  of  the  coefficient  matrices  for  the  fourteen  cases  of  Table  2. 
A  plot  of  the  of  the  minimum  separation  distance  divided  by  the  r  value  vs.  the  condition 
number  of  the  coefficient  matrix  is  shown  in  Figure  1.  There  are  15  points  for  each  of  the 
simple  knot  and  multiple  knot  algorithms.  The  first  fourteen  points  correspond  to  the  data 
in  Table  2,  in  the  order  given,  while  the  fifteenth  is  for  the  thinned  glacier  dataset,  GL, 
using  25  knots.  For  the  latter  data  the  underlying  function  is  unknown.  The  simple  knot 
algorithm  found  9  near-double  knots.  The  multiple  knot  algorithm  found  four  double  knots 
and  three  triple  knots. 

The  points  in  Figure  1  for  the  simple  knot  algorithm  have  been  labeled  with  their  number. 
Because  of  the  clustering  of  many  of  the  points  for  the  multiple  knot  algorithm,  only  those 
with  condition  number  great  than  105  have  been  labeled. 

What  can  be  observed  is  that  the  condition  numbers  for  the  multiple  knot  algorithm  are 
generally  one  to  two  orders  of  magnitude  smaller  than  those  for  the  simple  knot  algorithm. 
There  are  exceptions  to  this,  most  notably  the  datasets  3.1,  3.2,  3.7.  and  4.5,  where  the  be- 
havior is  reversed,  and  3.6  for  which  the  r  value  is  relatively  large  in  both  cases,  contributing 
to  the  size  of  the  condition  number. 

Finally,  for  one  of  the  parent  surfaces,  we  include  a  plot  of  the  underlying  surface  being 
approximated,  along  with  the  surfaces  obtained  from  the  simple  knot  algorithm  and  the 
multiple  knot  algorithm.  The  dataset  is  3.7,  the  surface  shown  in  Figure  2a  being  sampled 
at  100  points.  Figures  2b  and  2c,  respectively,  show  the  surface  reconstructed  using  the 
simple  knot  and  multiple  knot  algorithms  with  twelve  knot  points.  In  this  case  it  is  clear 
that  the  multiple  knot  algorithm  has  achieved  a  better  fit. 

4      Conclusion  and  Ideas  for  Further  Work 

We  have  formulated  an  algorithm  to  allow  coalescence  of  near-repeated  knots  into  multiple 
knots,  bringing  into  the  approximation  the  appropriate  derivatives  of  the  basis  function  in 
place  of  the  basis  functions  for  the  nearby  knots.  As  a  general  rule,  this  results  in  both 
better  approximations,  and  more  stable  computations,  although  individual  examples  do  not 
conform  to  this  ideal.  Our  computations  convince  us  that  the  occurrence  of  multiple  knots 
is  a  natural  phenomenon,  and  that  algorithms  for  the  problem  need  to  account  for  the 
possibility. 

The  tendency  toward  repeated  knots  leads  to  one  possible  idea  for  further  exploration. 
This  would  be  to  begin  initially  with  some  (or  all)  knots  being  double  knots  to  take  advantage 
of  the  apparently  greater  fitting  power  of  the  multiquadric  plus  its  directional  derivative.  We 
have  not  yet  had  the  opportunity  to  pursue  this  idea.  While  the  advantage  of  starting  with 
multiple  knots  seems  clear,  the  disadvantage  is  that  coalescence  of  these  knots  skips  odd- 
multiplicity  knots,  especially  triple  knots,  which  occur  regularly  in  our  computations. 

We  have  observed  the  tendency  toward  repeated  knots  with  the  LEASTSQ  algorithm 
and  another  (FMINS  -  see  [W085])  in  the  previous  paper  and  feel  the  repeated  knots  are 
not  occurring  as  an  artifact  of  the  optimization  process.  Nonetheless,  an  algorithm  that  took 


advantage  of  the  particular  problem  being  solved  could  be  more  efficient  and  possibly  avoid 
some  of  the  local  minima  being  found  by  the  general  minimization  routines. 
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Figure  Captions 


Figure  1:  Minimum  separation  distance  relative  to  r  value  vs.  condition  number  for  thirty 
problems.  The  o's  denote  simple  knot  cases,  numbered  from  1-15,  the  first  14  being  as 
ordered  in  Table  2.  Number  15  is  the  thinned  glacier  data,  GL,  with  25  knots.  The  x's 
denote  the  multiple  knot  cases  -  only  those  with  condition  number  greater  than  105  are 
labeled  due  to  space  considerations. 

Figure  2:  a)  The  parent  surface  from  which  the  data  was  obtained,  b)  The  surface  obtained 
with  12  knots  using  the  simple  knot  algorithm,  c)  The  surface  obtained  with  12  knots  using 
the  multiple  knot  algorithm. 
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